Orientational correlations and the effect of spatial gradients in the equilibrium steady 
state of hard rods in 2D : A study using deposition-evaporation kinetics 
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Deposition and evaporation of infinitely thin hard rods (needles) is studied in two dimensions 
using Monte Carlo simulations. The ratio of deposition to evaporation rates controls the equilibrium 
density of rods, and increasing it leads to an entropy-driven transition to a nematic phase in which 
both static and dynamical orientational correlation functions decay as power laws, with exponents 
varying continuously with deposition-evaporation rate ratio. Our results for the onset of the power- 
law phase agree with those for a conserved number of rods. At a coarse-grained level, the dynamics 
of the non-conserved angle field is described by the Edwards- Wilkinson equation. Predicted relations 
between the exponents of the quadrupolar and octupolar correlation functions are borne out by our 
numerical results. We explore the effects of spatial inhomogeneity in the deposition-evaporation 
ratio by simulations, entropy-based arguments and a study of the new terms introduced in the free 
energy. The primary effect is that needles tend to align along the local spatial gradient of the ratio. 
A uniform gradient thus induces a uniformly aligned state, as does a gradient which varies randomly 
in magnitude and sign, but acts only in one direction. Random variations of deposition-evaporation 
rates in both directions induce frustration, resulting in a state with glassy characteristics. 

PACS numbers: 64.60.Cn 
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I. INTRODUCTION 

An assembly of particles interacting via hard-core re- 
pulsion serves as a useful model for studying simple flu- 
ids, colloids, liquid crystals and many other soft matter 
systems. The analysis of such model systems helps in 
understanding the features of real systems, such as their 
phase behaviour, structural and dynamic properties. An 
important role is played by the anisotropy of shape of 
the constituent particles, which can range from thick 
elongated platelets to thin rods. Some examples of sys- 
tems in which the constituent particles show anisotropy 
are certain types of colloids, liquid crystals and protein 
molecules. In particular, rod-like particles are found in 
suspensions of the tobacco mosaic virus Q], nematic liq- 
uid crystals 2] and, recently, carbon nanotube gels 0. 
All these systems show very rich and characteristic phase 
behaviour. 

Rod-like particles have been modeled theoretically as 
ellipses 0, rectangles and spherocylinders [UIHOl with 
varying aspect ratios, a limiting case being infinitely thin 
hard rods or needles 8J. These systems exhibit a num- 
ber of interesting entropy-driven phase transitions which 
have been studied in two and three dimensions, usu- 
ally using simulations with number-conserving dynam- 
ics. On the other hand, there are a number of physical 
processes which involve adsorption (deposition) and des- 
orption (evaporation) of particles, which do not conserve 
particle number and which are important for some mono- 
layer growth processes. Adsorption and desorption are 
also important in the binding/unbinding of ligands to mi- 
crotubules, the interaction of proteins with DNA [tj [^} 
and many catalytic reactions. Finally, in recent exper- 
iments on assemblies of long objects (rice grains, thin 



metal rods) on a vibrating plate ^1 > individual particles 
jump off and return to the plate, leading ultimately to 
a state with interesting patterns. These considerations 
motivate us to study the deposition and evaporation of 
hard objects with rigid boundaries on a substrate. While 
a deposition-only system, of the type studied in random 
sequential adsorption 0, can end up in a non-evolving 
jammed configuration, with the addition of evaporation, 
the system eventually reaches an equilibrium steady state 
with a density governed by the rates of deposition and 
evaporation [12-17]. While most of these studies have fo- 
cussed on the kinetics of approach to steady state, in this 
paper, we are interested in the properties of the steady 
state itself. Specifically, we study the patterns formed 
due to deposition and evaporation of infinitely thin hard 
rods (needles) on a 2D substrate. Needles are a limiting 
case of rod like particles in the systems mentioned earlier. 
Though not directly applicable to any physical system, 
this is an important limiting case; the limit of zero width 
simplifies the problem by eliminating the aspect ratio as 
a parameter. The hard core constraint is enforced by re- 
jecting any deposition event which results in an overlap 
of needles. 

It is useful to recall some known facts about a sys- 
tem of hard needles with no externally imposed spatial 
inhomogeneities. This system shows a transition from 
a low-density orientationally disordered (isotropic) state 
to a high-density ordered state with nematic correlations. 
This transition, whose existence was pointed out by On- 
sager |18(, can be viewed as an outcome of the interplay 
between orientational and translational entropy of the 
needles; the ordered (aligned) state is preferred at high 
density since alignment leads to an increase of transla- 
tional entropy, albeit at the cost of orientational entropy. 
The nature of the orientational ordering is dimension- 
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dependent. In three dimensions, orientational long range 
order (LRO) sets in. A state with LRO would break the 
continuous symmetry of rotations, and is thus not ex- 
pected to occur in 2D, even though the Mermin- Wagner 
theorem cannot be generalized to this system [l^. In- 
deed, the simulation study of Frenkel and Eppenga 
on a system with a fixed number of needles confirms the 
absence of LRO in 2D, and finds a phase with power-law 
decays of orientational correlations, quite analogously to 
the XY model 20]. 

On a coarse-grained scale, the local orientation at loca- 
tion r and time t is specified by an angle field 9(r, t). The 
orientational correlation functions of interest are defined 



as 



gt{T,t) = (cos[t(fi(T,t) -0(0,0))]) 



(1) 



where £ is an even integer, and 9 and 9 + n represent the 
same state. Quadrupolar correlations are probed by £ = 
2, whereas higher values of £ correspond to higher mul- 
tipolarities. From numerical simulations, we find power 
law decays in both space and time: ge(r, 0) ~ r~ Ve and 
ge(0,t) ~ t~ 0e for both £ = 2 and 4. Our results for 
the static correlations conform to the Kosterlitz-Thouless 
theory for the onset of correlations, while our results for 
the dynamics show that their decay in time is governed by 
the Edwards- Wilkinson equation. We also study spatial 
variations of the deposition rate, and find strong effects 
on the nature of the ordering. We consider several types 
of variations: (i) a sharp change across a linear interface 
; (ii) a smooth linear gradient; (iii) a random variation of 
rates in one direction and (iv) random variation of rates 
in the plane. We find that the qualitative effect of spatial 
variations is to induce alignment of needles in the direc- 
tion of gradient, an effect with an entropic origin. In 
(i), the effect dies down slowly with increasing distance 
from the interface, but in cases (ii) and (iii) it results 
in a state with overall alignment along the direction of 
variation of the deposition rate. The random variation 
in (iv) induces frustration and the result is then a state 
with glassy features such as strong initial condition de- 
pendence and slow relaxation. 



II. 



MODEL AND PROCEDURE 



In our model, infinitely thin hard rods (needles) are 
added to a 2D substrate with area A with a constant at- 
tempt rate, and simultaneously removed randomly from 
the substrate with a specified rate. In a deposition at- 
tempt, the location of the centre of mass of the needle is 
chosen at random on the substrate, and the orientation 
angle is chosen at random as well. Let Td be the rate 
of attempted depositions per unit area per unit angle in- 
terval. An attempt is successful only if the depositing 
needle in question would not overlap with existing nee- 
dles on the substrate; otherwise it is rejected. 

During evaporation, a needle is chosen at random from 
those present on the substrate and then removed. If the 



total rate of such removals is R e we may associate a re- 
moval rate T e =-^j per unit area per unit angle interval. 
The ratio 



(2) 



of deposition to evaporation rates is the control param- 
eter in the problem. As we show below, k is related to 
the fugacity z — of an equilibrium grand canonical 
system. 

The model under consideration can be thought of as 
describing the adsorption and release of needle-like gas 
molecules on a substrate in contact with a gas reservoir 
with which it can exchange particles. The equilibrium 
state on the substrate is then described by the grand 
canonical (fiAT) ensemble, where /i , A and T are respec- 
tively the chemical potential, substrate area and temper- 
ature. Define scaled coordinates Sj = r^/L, where is 
the position of the i th particle on the substrate and L is 
the linear dimension of the system. The grand canonical 
partition function can be written plj as 



S = E Jn OVAT / <fai / ds 2 .... / ds N 

J ddi J d0 2 J d0 N e- pu (w>~»»fr>**,.-9x) (3) 

where A = ( 2 ! r? '_ ) 1 / 2 is the thermal wavelength which 
results from integrating over the momentum of needles, 
and U is the interaction energy for a configuration in 
which there are N needles with scaled centre of mass lo- 
cations (si, S2, ...sjy) and orientations {6%, 62, ....Ojsi). The 
corresponding equilibrium probability density_ of a con- 
figuration C = (si, s 2 , ...sjv, 9i, 9 2 , ....9n) is 



lty c 

m 



N 



(4) 



For our system of hard core needles, the interaction en- 
ergy U — > 00 when needles overlap, while [7 = when 
there is no overlap between any needles. Thus all allowed 
configurations with fixed N have equal weights. 

Deposition - evaporation dynamics involves changes of 
N. The evolution from a configuration C to a configura- 
tion C can be described by the master equation 



dP(C) 

at 



C) - P{C)W(C -> C')} 

p (5) 

The steady state of Eq. (JSJ, obtained by setting 9P q^ 
— 0, is in fact an equilibrium state if the condition of 
detailed balance 



W{C -> C) = P eq (C) 
W{C -> C) P eq {C>) 



(6) 



3 



is satisfied for every pair of configurations C and C that 
can be reached from each other. Now, let C denote an 
,/V-needle configuration and let C be the (N — l)-needle 
configuration obtained from C by removing a particular 
needle. Using Eq.(@} in Eq.©, we see that = j^. 
Thus, the steady state of deposition-evaporation dynam- 
ics is described by the grand canonical equilibrium state 
with 

where p = ^ is the areal density of needles. 

In our Monte Carlo studies we varied the control pa- 
rameter k in the range 1 to 40 and monitored the re- 
sulting density and orientational correlations. We used 
an L x L substrate (with L = 15 and 25) where L is in 
units of needle length. For each value of k we made mul- 
tiple runs, allowing up to 10 7 Monte Carlo time steps for 
equilibration. The Monte Carlo time t is defined as the 
number of attempts divided by L 2 . Averaging was done 
over 10 sets of independent runs and 1000 configurations 
from each run after equilibration. 

Figure [l (a)| sh ows the variations of the density with k, 
while Fig. |l(b)| sh ows p plotted against z/A 2 = pn. The 
inset in Fig. |l(a)| shows a marked change in the depen- 
dence of pi (k — 1) on k for k in the range 20-25. As we 
shall see, there is a transition to a phase with power law 
decay of orientational correlations beyond n = k c ~ 25, 
as illustrated by the representative configurations shown 
in Fig. [21 for different values of k. We turn to a quantita- 
tive study of orientational ordering in the next section. 




(a) 




200 300 400 

•> 

z/A"( = pK ) 
(b) 



500 



III. ORIENTATIONAL ORDERING 

A. Order parameter 

For a system of N hard rods in 2D, the nematic order 
parameter q is given by 



= ^(f> S (2^)> 



(8) 



where 9i is the angle made by the i rod with the ne- 
matic director, which itself has an orientation (f> with re- 
spect to a fixed reference X-axis. Both </> and q can be 
found by studying the tensor order parameter defined as 



1 N 

= ^(^Z[ 2u «^)"' 3 (*) 



(9) 



2=1 



where u a (i) is the a th component of u(i), the vector spec- 
ifying the orientation of the i th rod. The eigenvalues of 
Q a/ 3 are ±g, and the corresponding eigenvectors pick out 
directions along and perpendicular to the director orien- 
tation 4>. Insofar as there is no long range order in the 2D 



FIG. 1: (a) The variation of p (number of rods per unit area) 
with « shows a change of behaviour for k in the range 20-25. 
This is more prominently depicted in the inset which shows 
the variation of p/(n — 1). (b) Variation of p with z/A 2 . The 
inset shows the initial portion of the curve. 



needle system, q vanishes in the thermodynamic limit. In 
simulations on finite systems, though, q may appear to 
2(d) |, over short times. Tracking the 



be non-zero (Fig. 
onset of such an apparent value is not a reliable way to 
locate the transition point. 



B. Orientational cumulant of q 

A better indication of the transition point, and also the 
nature of the ordered phase, is provided by monitoring 
the probability distribution functions Pl{<1) of q, where q 
is the block-averaged value of the local order parameter 
with the system divided into blocks of finite size L [22| . 
A measure of the non-Gaussian character of Pl (q) is pro- 
vided by the reduced A th order Binder cumulant of q 



U L = 1 



(g 4 )L 

3<<Z 2 >! 



(10) 
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(d) « ~ 39 

FIG. 2: Snapshots of hard rod configurations at different val- 
ues of re. The substrate size is 15x15. Observe the formation 
of defects in the configurations (a), (b) and (c). 



where L is the linear size of the sub-system (block). Ul 
provides a useful diagnostic tool to monitor the ordering 
induced by varying re |24j . 

For L <C £, C/z, is expected to stay close to a fixed point 
value U*. Thus, the occurrence of a critical point with 
£ = oo can be identified by plotting [/ £ against re for 
various values of L, and looking for common intersection 
point. 

We analysed the Monte Carlo data of our model by 
monitoring Ul (Eq. llOjl . The analysis was performed 
as follows : We simulated a single large system of size 
K x K (K = 25) and divided it into subsystems of size 
L x L, thereby having total M 2 number of subsystems 
with M = K/L [2i|] . Then, M was incremented in inte- 
ger steps starting from 1 and Ul was estimated for those 
subsystem sizes L where a good analysis is possible. Con- 
sequently, we did not consider very small or very large 
values of M. Also, the curve for M = 12 lies anomalously 
low and was not included. The number of subsystems 
(M) we use for estimating the cumulant range from 8 to 
20. Fig. shows the variation of Ul as a function of re 
for various values of L. 




K 

FIG. 3: Orientational cumulants Ul as a function of 
deposition-evaporation ratio re for various subsystem sizes L 
= K/M. Note the collapse of the curves beyond re c ~ 25, 
pointing to the occurrence of a power law phase. 

Below re c ~ 25.8, the curves are separate and distinct, 
but they collapse at re c , indicating the onset of ordering 
(Fig. [3J . Moreover, the curves seem to stay collapsed for 
re > re c suggesting that £ remains infinite in this phase, 
i.e. this is a phase with power law decay of correlations. 
Corroboration of this is provided by directly monitoring 
the correlation functions as described below. 



C. Orientational correlation functions 

Let us define a general orientational correlation func- 
tion ge(r, t) = {cos[£(9(r, t) — 8(0, 0))]) where I is an even 
integer. We studied static and dynamical properties by 
investigating g t (r,0) and g^(0,t). We calculated spatial 
correlations by forming circular bins around each rod in 
turn, computing gi(r, 0) for each bin, repeating this pro- 
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cess for all rods in the configuration and averaging over 
all rods (see Fig. 0J. The dynamical correlation func- 
tion, ge(0,t) was calculated by coarse-graining. The sys- 
tem was divided into a number of small cells (lxl) and 
an average value of orientation was assigned to each cell 
by averaging over the orientations of those needles whose 
centres of mass lie in the cell. The value of ge(0,t) was 
computed using this average value over each time frame, 
and averaging over all the c ells (s ee Fig. . The initial 
drops of the curves in Figs. 



4(b) and |5(b)j are sensitive 



to the size of the cell used, while the power laws seen at 
the larger times do not depend on the cell size. 

In the nematic-like phase, i.e. for k beyond k c , 
the correlations decay algebraically gi(r, 0) ~ r~ ve and 
ge(0,i) ~ t~P e . There are pronounced finite size effects 
which lead to a flattening of the curves for r as L/2, 
limiting the range over which the power law behaviour 
extends. We found that the values of the exponents rje 
and [3i vary continuously with k as shown in Table I. 
The estimations were done over 10 independent sets of 
configurations, each averaged over 1000 configurations. 
For k ~ k c , we observed that 772 — 0.23 ± 0.03 close to 
the predicted Kosterlitz-Thouless value 0.25. Our results 
for the static case, agree with those reported by Frenkel 
and Eppenga for the case of a fixed number of hard 
rods on 2D plane. We confirm that at the critical point 
k = k c , the mean density is ~ 7 

For k > k c , it was observed that exponents ob- 
tained from static correlations gi (r) and gi(r) are related 
through j]2 — ?74/4. It was also found that the exponents 
derived from the temporal correlations (72 (t) and g^ (t) 
are related in a similar way, i.e. /?2 — $1/4. Further, 
the ratios r\ilfii — 2.0 for I = 2 and 4, implying that the 
dynamical exponent Zd yn is 2. These observations can be 
understood on the basis of the simple model discussed 
below. 

In order to model the dynamics we note that the 
stochastic evaporation and deposition events change the 
local value of the coarse-grained angle field 8 in a noisy, 
diffusive way. In the discussion below, we take the an- 
gle to be an unconstrained variable running from -00 to 
+00 with (6 + ri7r) denoting the same needle orientation 
for integer n. We consider a simple phenomenological 
equation 



86 
dt 



K\7 2 8 + £ 



where £ denotes white noise which satisfies 



(11) 



(12) 



where B is a constant. This is of the same form as 
the Edwards- Wilkinson equation |2^|, which describes 
the evolution of a fluctuating interface. In our context, 
Eq. (|ll|l follows from the symmetric form of the Frank free 
energy F = h K j (\78) 2 d8 on using the phenomenologi- 
cal Langevin equation ^ = + £. A more complete 
description would involve coupled equations for the non- 
conserved density and orientational fields. From Eq. i|ll|) 



and i|12fl it follows that 
{Q(r + r',t + t')6[r,t)) 

Setting t' = we find that 
{[8{v + v',t)-8{r,t)] 2 ) 



B f d 2 k 



64ir 5 K / k 2 



.gik-r' e -Kk 2 t' 



(13) 



B 



32tt 5 A' 



x27iTn(r) (14) 



which, using the Gaussian property of 8, further implies 

(cos{£(8{r + r', t) - 8(r, t))}> ~ (15) 

where rje — £ 2 B/32ir 4 K. The measured values of rjg can 
be used to find K/B, whose value is included in Table I. 

Similarly, setting r' — in Ea. lfH^ we find the auto- 
correlation function 

(cos{£(8(r, t + t')- 8(r, t))}> ~ t~ 01 (16) 

where — £ 2 B / 64tt 4 K . Thus, for all k > k c the ratio 
of i]4 to 1J2 (and /J4 to f3 2 ) is expected to be 4; as we have 
seen above, our numerical results confirm this. Also, we 
find the dynamical exponent Zd yn — We/fli is — 2.0 . 

IV. EFFECT OF INHOMOGENEOUS k 

In this section, we explore the effects of having a spa- 
tial variation of the deposition-evaporation rate ratio k. 
In a physical system, such a variation could arise from 
the variation of the chemical potential or substrate tem- 
perature from one spot to another as their local values 
could influence the local detachment rate, as seen from 
Eq. J7J). As expected, such changes in k induce a spatial 
variation of the density; more interestingly, they have a 
strong effect on the orientational order as well. We ex- 
plore these effects by considering several types of spatial 
variations of k. 

In parallel to the discussion in section II, let F^ic) 
and r e (x) denote the deposition and evaporation rates at 
point x in the plane, and let k(x) = Let C be the 

configuration reached from configuration C by removing 
the rod at x m and let P(C) and P(C) be the weights 
of the respective configurations in steady state. We can 
check that the condition of detailed balance is valid when 
P{C) has the product form P(C) = Hiz(xi) where z{xi) 
A 2 p(xi)n(xi) is the local fugacity at the location of the 



;lh 



rod. If C is obtained from C by evaporating the m 



Hi 



rod, then evidently P(C) — P{C) x z(x m ). Recalling 
that W(C -> G) = T d {x m ) and W(C C) = T e (x m ), 
we see that the condition of detailed balance is valid. 

Thus the system reaches an equilibrium steady state 
which, however, is inhomogeneous in density, due to the 
non-uniform position dependence of k. The nature of 
orientational order depends strongly on the the way in 
which k is specified to vary over the plane. Below we 
consider several types of variations 

(i) A single interface separating low and high k regions, 
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(a) 



(b) 



FIG. 4: Log-log plots of (a) g 2 (r,0) = (cos[2(8(r,t) - 0(0, t))]) and (b) g 2 (0,t) = (cos[2(0(r,i) - 0(r,O))]) showing the static 
and temporal behaviour respectively of <?2(r, i). Data is shown for systems of sizes 15x15 and 25x25, and different values of k. 
The distance r is in units of rod length and time t is in Monte Carlo time steps. 
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(a) 



(b) 



FIG. 5: Log-log plots of (a) ff 4 (r,0) = (cos[4(0(r,t) - 0(0, t)))) and (b) y 4 (0,i) = (cos[4(0(r, t) - 0(r,O))]> showing the static 
and temporal behaviour respectively of g^{r,t), for the same system sizes and values of k as in Fig. 0] 



k{x) = «i, for £ <4, 
= «2, for x >-j. 

(ii) A uniform gradient in re across the substrate, 

k(x) = ki (1 + a f ) 

(iii) Random variation of ft in the X-direction only, 

k(x) — K% + 5k(x) 

where 8k(x) < k% is a random function of a;. 

(iv) A random binary assignment of k on a grid on the 
2D substrate 

In the first three cases, periodic boundary conditions 
were applied in Y-direction and open boundary condi- 
tions along X. In the last case, open boundaries were 
used in both directions. In all the cases considered, the 



ranges of k values were chosen to be above k c , the critical 
value in the uniform case. Our findings are as follows : 



A. Ki — K2 interface 

Here, a uniform value k± operates upto half-way across 
the 2D plane along the X-direction, while k — k,2(> Ki) 
in the remaining half. In the vicinity of the interfacial 
region, the rods are observed to orient in the direction 
the of k gradient i.e. perpendicular to the interface (see 
Fig, [gag. 

This can be understood on entropic grounds. That 
arrangement of rods is favoured which maximizes the en- 
tropy. By symmetry, the preferred average orientation 
of rods should be either (a) parallel or (b) perpendicu- 
lar to the interface. Consider those rods in the high-K 
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TABLE I: k dependence of exponents rjt and pi for I = 2 and 4 . The estimated error is indicated in the brackets. 



K 


P 


r\2 


V* 




At 


K/B 


26.0 


7.3 


0.23 (0.03) 


0.87 (0.06) 


0.11 (0.04) 


0.31 (0.06) 


5.3 xl0" a 


32.3 


9.4 


0.098 (0.006) 


0.41 (0.02) 


0.037 (0.007) 


0.16 (0.04) 


12.8 xlO" 3 


39.0 


11.4 


0.059 (0.003) 


0.25 (0.01) 


0.022 (0.003) 


0.09 (0.01) 


21.4 xl0~ 3 



half, whose centres lie very close to the interface (within 
half a rod length) so that part of such rods can reach 
into the low density side. Small variations in the angle 
of each rod would contribute to the entropy, but these 
are limited by the presence of other rods. A horizontal 
average alignment allows the rods to sample a less dense 
environment, and thus be subject to fewer constraints, 
on one side. Thus, option (b) would be preferred over 
(a) . The effect of interface- induced horizontal alignment 
is felt for some distance awa y from the interface on both 
sides. This is evident in Fig 6(a) which shows a steady 
state configuration in a system of size 25 x 15 (in units 
of rod length), with n\ = 30 and K2 = 50. However, for 
a large enough size, the system reverts to a power-law 
phase in the region far from the interface as observed in 
the uniform k case. The correlation function decays as a 
power law in the bulk, away from the interface, as shown 



in Fig 6(b) 



We also considered the case with two values of k, using 
periodic boundary conditions in both directions which 
leads to two interfaces (Fig. 0). The figure shows a 80 
x 15-sized system with periodically linked left and right 
quarters with n\ = 30, while the middle half has K2 = 50. 
Evidently, this system too shows interface-induced hori- 
zontal alignment. This geometry admits of an interesting 
limit where the entropy driven alignment is particularly 
clear. On shrinking the width of the central region to 
zero, at the same time taking the limit n\ —> 0, we obtain 
a ID model as a limiting case. In this model, deposition 
and evaporation moves are allowed with needle centres 
constrained to lie on the line. Needles are found to ori- 
ent preferentially perpendicular to the line (see the inset 
in Fig. |8J). The reason is evident. If the mean orientation 
of the director is perpendicular to the substrate, needles 
have the largest leeway to make angular excursions about 
the mean - i.e. the rotational entropy is then the largest. 
The variation with k of the density and order parameter 
q = {cos2(j)) where <p is the angle made with the direction 
perpendicular to the line, is shown in Fig. [S] 

A similar effect should also lead to rods getting aligned 
horizontally if they are close to an open boundary in the 
2D system. That this is so can be seen in Fig. [§1 which 
shows a system with uniform k = 27 and open boundary 
conditions along the X-direction. 




(a) 




FIG. 6: (a) Snapshot of a typical hard rod configuration with 
a single n\ — K2 interface. The system size is 25 x 15 and ki = 
30 (left half) and K2 = 50 (right half). Boundary conditions : 
Open (along X) , Periodic (along Y). Notice the difference in 
density in the two halves. (b)This plot shows the decay of the 
orientational correlation function (72 (y) = (cos[2(#(y) — #(0))]) 
calculated for a pair of points in the same vertical strip of unit 
width. Curves from bottom to top correspond to different 
strips in two halves in the configuration (a). The inset shows 
52(1") = (cos[2(6(r) — 0(0))]) measured radially for a box of 
size 6x6 which is positioned at the centre of each of the left 
and right halves of the same configuration. The distances r 
and y are in units of rod length. 
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FIG. 7: A dual-interface configuration in an 80 x 15-sized 
system with periodic boundary conditions in both directions. 
The middle half of size 40 x 15 and with K2 = 50 separates 
two quarters, linked at the boundary, each of size 20 x 15 and 
having Ki = 30. 



FIG. 9: A typical configuration of size 25 x 15 with uniform 
k — 27 and open boundary conditions along the X-direction. 
Free boundaries induce alignment which propagates some dis- 
tance into the bulk. 
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FIG. 10: A typical configuration for a 25 x 15-sized system 
with a uniform n gradient, with kl = 32 and hr = 50 at the 



two ends respectively. The horizontal alignment induced by 
the gradient is evident. 



FIG. 8: Variation of the density(p) and order parameter^) 
with k for the ID model which a limiting case of the 2-D 
Ki — n-2 model for Fig. |7| As shown schematically in the 
inset, the preferred orientation of needles is perpendicular to 
the line. 



B. Uniform k gradient 

This case is related to the discussion above, as the lin- 
ear increase in k may be viewed as a continuous succes- 
sion of interfaces from one end to the other. Since each 
interface induces an alignment of rods across it, this re- 
sults in overall alignment of the rods in the system (see 
Fig. I10[l . Note that the alignment is not an outcome of 
spontaneous breaking of orientational symmetry, as the 
gradient in n singles out a direction in space. We checked 
that the horizontal alignment is not tied to the aspect ra- 
tio of the container by simulating a system size 10 x 40, 
and observing overall horizontal alignment of rods in the 
steady state. 

Figure ITU1 shows a steady state configuration in a sys- 
tem of size 25 x 15, with n varying linearly from a 
value kl — 32 at the left end to a value kr = 50 at 
the right end. In our simulations, the system was equi- 
librated for 10 7 MCS and 10 4 post-equilibration config- 
urations were used to calculate averages. We studied 



the spatial and dynamical behaviour of the orientational 
correlation function g2(r,t) — (cos[2(9(r,t) — 9(0,0))]). 
Since the system is inhomogeneous along the X-direction, 
the substrate was divided into vertical strips, each hav- 
ing width of a rod length, and each strip was studied 
separately. The Y-density of needles inside each strip 
was uniform though the density varies from strip to 
strip. We monitored the correlation function 32(2/) = 
(cos[2(9(y + y ) - 9(y ))}) (see Fig. [TTJ, and found that 
52(2/) decays exponentially to a non- vanishing constant 
value qo which differs from strip to strip. For the sys- 
tem under study, <?q varies in the range 0.68-0.77 over 
the strips. 

The dynamical correlation function, i.e. 92(t) = 
(cos[2(9(t) — 9(0))]) was calculated by coarse-graining. 
The system was divided into number of small cells (2x2) 
and an average value of orientation was assigned to each 
cell by averaging over needles in it. The plot of gi (t) is 
shown in Fig. for cells at different values of X. Each curve 
shows an exponential approach to a non-zero constant 
value. The behaviour of both the spatial and dynamical 
orientational parts of the correlation functions indicates 
a phase with overall orientational alignment. 
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FIG. 11: Correlation function in a 25 x 15-sized system with 
a uniform re gradient. Shown is the log-normal plot 32(3/) = 
(cos[2(6(y) — 6(0))])— ql for pairs of points in the same vertical 
strip of unit width. Curves from bottom to top correspond 
to different vertical strips in the order of increasing re. It is 
evident from the plots that g~2(y) exhibits exponential decay. 
The distance y is in units of rod length. 



Figure ITSI shows a steady state configuration obtained 
with a quenched random variation of k(x). The system 
was simulated by varying re randomly around a value of 
k such that k ± 5k(x) > re c , where Sn(x) denotes ran- 
dom variations along X-direction. We used K ~ 32 and 
5k = 2.0. The correlation functions 172(2/) and 172 (t) be- 
have similarly to the uniform gradient case (see Fig. 1141 
andlSI). Thus, this case also yields a phase with overall 
orientational alignment. 




FIG. 13: A typical configuration for a 25 x 15-sized system 
with a re gradient achieved by varying the value of re randomly 
around 32 in the X-direction only. The horizontal alignment 
induced by the gradient is evident. 
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FIG. 12: Dynamical correlation function in the system of Fig. 
[TDlandini Log-normal plots of g 2 (t) = (cos[2(6(t) - 9(0))]) - 
go as shown, with curves from bottom to top correspond to 
different cells with increasing « (refer to the main text for 
details). The time t is in Monte Carlo time steps. 



C. Random variation of k(x) 

We have seen that a uniform gradient in k results in 
an orientationally ordered state. However, the argument 
for ordering does not depend on the gradient being con- 
stant in magnitude or sign. Thus, if n varies randomly 
(with k > k c ) along the X-direction but is uniform along 
Y, the resulting state once again should exhibit horizon- 
tal alignment with needles aligned along the X-direction. 
The resulting state can once again be viewed as contin- 
uous succession of interfaces and should display overall 
alignment. 




FIG. 14: Correlation function in a 25 x 15-sized system with 
a random re gradient in the X-direction only. Shown is the 
log-normal plot 32(2/) = (cos[2(8(y) — 9(0))]} — qo for pairs of 
points in the same vertical strip of unit width. Curves from 
bottom to top correspond to different vertical strips in the 
order of increasing re. It is evident from the plots that g2(y) 
exhibits exponential decay. The distance y is in units of rod 
length. 



D. 2D random binary distribution of re values 

In this case, the fugacity is set inhomogcneously in 
a quenched disordered fashion, so that the tendency to 
align locally along gradients results in competing pat- 
terns of order, i.e the system is frustrated. The resulting 
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FIG. 15: Dynamical correlation function in the system of Fig. 
OandfTl. Log-normal plots of g 2 (t) = (cos[2(6(t) - 6(0))}) - 
§q as shown, with curves from bottom to top correspond to 
different cells with increasing k. The time t is in Monte Carlo 
time steps. 



state has glassy features and contains domains of differ- 
ent orientations (see Fig. Ilfijl [2f|. 

In our simulations, we divided the substrate of size 25 
x 25 into a grid of unit length squares. Each square 
was randomly assigned a k value either 27 or 50 (both 
greater than n c ). The random k gradient across square 
edges generates local disorder, which can disrupt the ori- 
entational order and result in destruction of orientational 
alignment on the scale of the system size. The effect of 
quenched random disorder due to orientational random- 
ness of cross-links in a system of nematic elastomers has 
been studied earlier (2tJ and the model was reported to 
have spin-glass like behaviour. In our model, the disorder 
emerges from the randomness in the spatial distribution 
of k values. We find that the spatial correlation (r) 
decays exponentially to zero (Fig. I17fl whereas the dy- 
namical part giit) seems to decay in an algebraic manner 
to a non-zero value (Fig lTH)) . 

The behavior is suggestive of a glassy system which 
is disordered in space but relaxes slowly in time. More- 
over, it was also found that with same quenched disorder 
arrangement, different initial conditions lead to different 
states. Typical configurations in each of these states are 
shown in Fig. ^3 which is a glass-like feature. 



V. CONTINUUM DESCRIPTION 

As discussed above, our simulations show that a spa- 
tially inhomogeneous deposition-evaporation ratio k can 
induce nematic order and other interesting orientational 
patterns in the equilibrium state of a system of hard nee- 
dles. It is interesting to ask whether these effects can be 
captured within a phenomenological coarse-grained de- 
scription based on including symmetry-allowed terms in 
the free energy. In the context of liquid crystals, such an 
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FIG. 16: Snapshots of hard rod configurations for random 
binary distribution of k values 27 and 50 on the substrate of 
size 25 x 25. Representative configurations characterizing two 
different states ((a) and (b)) which are reached from different 
initial conditions, for the same k distribution. 



approach has proved successful in studying large-distance 
phenomena, including the effects of walls and other inho- 
mogeneities Q- Below we sketch such a description for 
our system of interest j2^. Besides showing that gradi- 
ents in the deposition-evaporation rates lead to orienta- 
tional ordering, the treatment suggests the occurrence of 
local splay. 

Let us define a nematic director field n(r) to describe 
the local coarse-grained value of the orientation of nee- 
dles (evidently, h(r) and —h(r) describe the same config- 
uration) . In the absence of externally imposed inhomo- 
geneities, spatial variations of n(r) lead to a free energy 
described by the Frank form 



Fk = J d 2 r[^(V •n) 2 + |(nx(Vx n)f] (17) 

The two terms describe, respectively, contributions of 
splay and bend to the free energy; there is no contri- 
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by Vk in the terms in Ea. (|17|) to get 




FIG. 17: Evidence for exponential decay of spatial correla- 
tions in a system with random binary distribution of k. The 
figure shows log-normal plots of 32(1") = (cos[2(6(r) — 0(0))]). 
The parameters are the same as in Fig. 1161 and the curves 
correspond to the different steady states evolved from two 
different initial conditions. The distance r is in units of rod 
length. 




FIG. 18: Evidence for a power law decay of temporal correla- 
tions in a system with random binary distribution of n. The 
figure shows log-log plots of c/2(t) = (cos[2(8(t) — #(0))]) — go- 
The parameters are the same as in Fig. Il(il and the curves 
are obtained in the different steady states reached from two 
different initial conditions. The time t is in Monte Carlo time 
steps. 



Fj = J d 2 ry [(Vk • n) 2 ] + J d 2 r'^-[(n x (Vk x h)) 2 ] 

.(18) 

In addition to these terms, which are quadratic in 
Vk, one can also construct scalars which involve Vk lin- 
early 29] 

F L = J d 2 r^-Vn-[h(V -n)} + J d 2 r^-V n-[(n ■ V)n]. 

(19) 

Symmetry considerations alone do not suffice to deter- 
mine the values of the coefficients K\, K 3 , J\, J3, L\ and 
L3. Their density dependences can be found on noting 
that a change in k induces a change in density, thereby 
influencing the elastic energy. To incorporate this, we re- 
place h in Eq. Q17[l by (p(r)n/ po) where p(r) and po are 
the local and average densities respectively. Appendix A 
contains the resulting form of the free energy F p and the 
values of the coupling constants. 

Let us turn to the consequences of the new terms. With 
a uniform spatial gradient, Vk = ax (case (ii) above), Fj 
induces overall alignment of needles. To see this, consider 
Fk+Fj. Evidently, Fk is minimized by any arrangement 
in which h is uniform in space while Fj is orientation 
dependent. Writing n — xcosip + ysiruf), we find 



Fj = a 2 (J\COS 2 4> + J 3 sin 2 



(20) 



F K 4 
with 



Fj is minimized by having an aligned state, either 
i> = (if J 3 > Ji) or with cp = § (if J 3 < J x ). Our 
numerical results, supported by the entropic considera- 
tions given above, show alignment in the direction of the 
gradient, implying A J = J3 — J\ is positive. 

Now consider the effect of Fl. Setting Vk = ax in Eq. 
(fTU|l . we find 

dcf) 

F L = —a(Li + L 3 ) sincj) cos<j> — 

ox 

dd) 

+ a(Licos 2 (f> — L 3 sin 2 (j)) — (21) 

oy 

(Fk + Fl ) can be minimized on noting that each of the 
Li and L3 terms is an eigenfunction of the Frank elastic 
matrix. The result is 



5^ -Li 2 L 3 . 2 
— = —r^cos <p + —Bin <p 
oy K 3 



(22) 



bution from twist in 2-d. 

Inhomogeneities in deposition-evaporation rates lead 
to spatial gradients Vk, which imply new terms in the 
free energy. These terms consist of scalars involving Vk 
and n, respecting invariance under h <-> — n. Two such 
scalars are obtained by replacing the gradient operator 



d~x = { K- l + K; ) C0S(t> sm<t) (23) 

These terms describe a spiralling tendency of the director 
in space. 

The full problem involves minimizing Fj + Fk + Fj, . 
If Fj is dominant, the director is primarily aligned along 
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the gradient implying <f> is small. Eas. l|22|) and ill'MI) then 
reduce to d<j)/dy « —L\/K\ which describes a spir ailing- 
director; further, Fj restricts angular excursions to be 
at most <fio oc -^= . Thus the predicted state is one with 
overall alignment along the gradient, but with local splay 
structures, each with a small opening angle ~ 2<pQ. This 
picture is borne out by our simulations. 

For case (iii) in which k(x) varies randomly, we see that 
Fj = a 2 (Jicos 2 (f) + J^sin 2 (j)) where a 2 is the spatial av- 
erage of the mean squared gradient. As for case (ii), the 
free energy is minimized by having a state with align- 
ment along the gradient, as observed in our simulations, 
provided J3 > J\. 

In case (iv), the gradients that appear in Eq. (jlSjl are 
random in direction leading to frustration in the arrange- 
ment of needles. Equations (|1 7|) . H18|) and H19|) provide a 
starting point for a theoretical description of the glassy 
state that results. 



VI. CONCLUSION 

In summary, we have studied orientational ordering in 
a 2D grand canonical system of hard rods using deposi- 
tion and evaporation moves. The control parameter is 
the ratio n of deposition and evaporation rates, which 
controls the density. The system with uniform k displays 
a transition from an isotropic phase (for k < k c ) to a 
phase characterized by algebraically decaying static and 
dynamical orientational correlations for k > k c . Further, 
the values of the critical exponents and the behaviour of 
the orientational cumulant are consistent with Kosterlitz- 
Thouless theory. The numerical results for the dynamical 
correlation functions are described by a phenomenolog- 
ical Edwards- Wilkinson equation for the non-conserved 
orientational field. 

Our principal results pertain to the new behavior in- 
duced by having a position-dependent n, and hence a 
space-varying density of rods. An anisotropic variation 
of k (say along the X-direction only) results in needles 
aligning along the k gradient. This was understood by 
first considering the effect of an interface separating re- 
gions with two values of k. Entropic considerations lead 
the needles to align normal to the interface, i.e. along the 
gradient. From another point of view, K-gradients lead to 
new terms in the Frank-like free energy and these in turn 
imply orientational ordering. Finally in a system with 
quenched disorder corresponding to spatially random k, 
we found indications of orientationally frozen states with 
glass-like characteristics. It would be useful to have a 
better characterization and understanding of this glassy 
state. 

The mechanism behind gradient-induced orientational 
ordering is simple: spatial variations of k induce varia- 
tions in needle density; and an average alignment of nee- 
dles along the gradient is preferred as this leads to an en- 
hanced entropy of rotational excursions around the mean. 
In effect, the K-gradient thus behaves like an external 



field acting to produce nematic order, ultimately due to 
the strong coupling between spatial and orientational de- 
grees of freedom in the needle system. Gradient-induced 
ordering effects should be present in three-dimensional 
systems as well. In 3-D, the Onsager mechanism for ne- 
matic long range order would predict ordering for values 
of a uniform k exceeding a critical value n c . The addition 
of a uniform K-gradient would be expected to lead to a 
nonzero value of nematic ordering for all values of k, and 
to enhance its value for k > k c . It would be interesting 
to test this prediction, and have a quantitative measure 
of gradient-induced ordering in 3D. 
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APPENDIX A 

To incorporate the effect of spatial variation of the den- 
sity, we write (p(r)n/po) in place of the director n in the 
Frank free energy Fk of Eq. I|17|) . The resulting expres- 
sion F p for the free energy can be written as 

F P = I d 2 r[^(V • pnf + ^(pn x (V x pn) f] (24) 
J z Po z Po 

The expansion of the integrand involves ten terms : 
[& p 2 sin 2 <b + (J% p 2 ) p 2 cos 2 ^} {d^/dxf 

+ 

+ ° 

[gi cos 2 <\> + (fi p 2 ) sin 2 cj)} (dp/dx) 2 

+ 

[|| ^ + ( J| p2) cos 2 0] {dp/dy f 
+ ° 

+ (#§ P 2 ) } 2p sine/) cos<p (d(j)/dx) (dp/dx) 

+ 

[|| - (|| p 2 ) ] 2p smtp cos<t> (d<j>/dy) (dp/dy) 
+ ° 

[#| - (if p 2 ) 1 2 sm( t> cos ^ ( d p/ dx ) ( d Pl d y) 
+ 

+ (#i P 2 ) } 2p 2 sin(j> cos(j) (d<p/dx) (d(f>/dy) 

+ 

[=&■ sin 2 (j) - (fi p 2 ) cos 2 <j) ] 2p m/dx) (dp/dy) 

+ 

cos 2 ^ + (|| p 2 ) sin 2 ^ ] 2p (dp/dx) [d<t>/dy) 
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In case (ii) where k varies linearly with x, the den- 
sity gradient is non-zero only along the X-direction, and 
vanishes along Y, so the terms involving dp/dy do not 
contribute. We can then read off the density dependence 
induced in the elastic constants in terms of the original 
Frank's constants as 



K[{P) = P 2 and K' 3 (p) = p- 



_ Ka 



Now, comparing the third term of the expression with 
Eq. H18|l and writing dp/dx = a((p) where £(p) = dp/ 8k, 
we have 



Similarly, grouping the fifth and the tenth terms together 
and comparing with Eq. 1)19(1 , we obtain 



L x (p) = ^2 P C(p) 
Po 



L 3 (p) = (-^p 2 )2p((p). (26) 
Po 



Mp) = ^(C(p)) 2 

Po 



Up) = ^p 2 (C(p)) 2 (25) 

Po 
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